It was a sunny day in Istanbul. The chief marketing officer (CMO) of Silverline, an award-winning Turkish supplier of built-in appliances for cooking, cooling, and dishwashing, was on a video call with the sales manager (SM) from his office. In a worried voice, the SM said:
“Look – we have great appliances. All of our retailers say that customers love them. But our favourite store in Istanbul complains about the same issue over and over again. They just can’t get enough of our products – in particular, the extractor fans. So they end up pointing customers to competitors’ products. This is a big opportunity loss for us. We should do something about this.”
A silence fell upon the room and the CMO said:
“You are right. My heart bleeds when I see our great products fail to take their deserved share in the market. We should keep good ties with all the retailers as they are the gatekeepers to the market. But perhaps we also need an analytical demand forecasting solution at store level. What do you think?”
The SM responded excitedly:
“Your point on having an analytical solution is great! As far as I know, the use of analytical approaches in our sector is not very common. If we have an automated demand forecasting tool and apply it for each store, we may have a huge advantage compared to our competitors.For stores with the stockout issue, we can try to predict the future demand in a similar store and use that as an estimate”
The CMO was happy that her voice was echoed by the SM:
“What a great idea! As far as I know, the use of analytical approaches in our sector is not very common. If we had an automated demand forecasting tool and applied it for each store, it might well give us a huge advantage over our competitors. For stores with the stockout issue, we can try to predict future demand in a similar store in terms of demographics and other characteristics, and use that as an estimate.”
The CMO was pleased that her suggestion resonated with the SM: “Perfect! Why don’t we dive right into this? I will talk to the production team and try to find out how to build a data analytics team that can provide some data-driven solutions.”
The following day, the CMO had a meeting with the production team and towards the end of the meeting said: “I am so excited! Let’s build this analytics team together. The team can use some forecasting tools and tell us how demand will look like in the next several months. This way, we can support the stores that have stock issues in relation to our extractor fans.”
Not long afterwards, the analytics team started to work with a dataset from a similar store to the one highlighted. The data covered the historical monthly sales of an extractor fan over a time span of 122 months.
The question was straightforward:
What would the future demand be for this product over the next year?
To be able to address the question in the case study, we will use three time series modelling approaches:
Autoregressive Integrated Moving Average (ARIMA)
Simple Exponential Smoothing (SES)
Holt-Winters' Model
Importing the Python libraries¶
Here, we will import the essential libraries required for our case study.
# Required Python libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import seaborn as sns
The following code chunk sets up the working directory.
#Find the current working directory
import os
cwd = os.getcwd()
print(cwd)
#Change the working directory
os.chdir('C:/Users/gyildiri/Downloads/chapter 9')
print(os.getcwd())
C:\Users\gyildiri\Downloads\chapter 9 C:\Users\gyildiri\Downloads\chapter 9
Data visualization¶
We may want to have a quick look at how the data set and variables look like. As an example, by running the code chunk below, we can generate a table that provides a brief preview of the first 10 rows of our dataset.
data = pd.read_excel('silverline_data.xlsx')
# Displaying the first 10 rows and 2 columns of the data
data_view = data.iloc[:10, :2]
print(data_view)
Time Sales 0 1 116 1 2 184 2 3 211 3 4 228 4 5 283 5 6 219 6 7 178 7 8 141 8 9 175 9 10 146
To get a feel of the data patterns, it is usually a good practice to perform a visual inspection. Here, we will obtain the simple time series plot of the raw data using the plt.plot function.
plt.plot(data['Sales'], color='blue')
plt.title('Sales')
plt.show()
Question: What is your interpretation of this plot? Is there a seasonal or trend pattern in the data?
We are ready to proceed to model estimation and prediction. We will perform the following steps for each of the models:
Split the data into training and test set, and estimate a demand forecasting model using the training set;
Predict sales in the test set and evaluate forecast accuracy via the root mean squared error (RMSE).
Re-estimate the forecasting model using the whole data set and generate sales predictions for the next 12 months.
Finally, we will compare the models based on the RMSE metric and pick the best model.
# Log transformation of Sales data
data['LSales'] = np.log(data['Sales'])
# If your variable contains zero values, you should add a small number (e.g. 1) before you do the log transformation:
#data['LSales'] = np.log(data['Sales'] + 1)
Note that our sales data do not contain zero values. However, in some datasets we may observe variables data that do take zero values. In such cases, we need to add a small number (e.g. 1) to each observation in the data before we perform the log transformation. In doing so, we avoid having log(0), which is undefined.
Next, to examine the SES model’s predictive ability, we need to split the data into training (in-sample) and test (out-of-sample) sets. To do this, we apply the 80/20 rule, using the first 80% of the observations for the training set and the remaining 20% for the test set. Given that we have 121 observations in total, we should use the first 96 observations for the training set, and the remaining 25 observations for the test set. The following code chunk splits our data into training and test sets.
# Simple Exponential Smoothing
ses_train = data['LSales'][:96]
ses_test = data['LSales'][96:]
Parameter tuning and estimation¶
To find the best alpha (smoothing parameter), we can ask Python to take different values between 0.1 and 0.9 in steps of 0.01, estimate the SES model, and calculate the corresponding RMSE metric for each alpha. Then we can plot all the RMSE values against values of alpha. Finally, we can choose the alpha parameter that yields the lowest RMSE. The following code performs parameter tuning for the SES model and generates the output below.
from statsmodels.tsa.holtwinters import SimpleExpSmoothing
from sklearn.metrics import mean_squared_error
# Parameter tuning to find the optimal alpha
alphas = np.arange(0.1, 0.91, 0.01)
ses_rmse_values = []
for alpha in alphas:
model = SimpleExpSmoothing(ses_train).fit(smoothing_level=alpha, optimized=False)
ses_predictions = model.forecast(len(ses_test))
ses_rmse = np.sqrt(mean_squared_error(ses_test, ses_predictions))
ses_rmse_values.append(ses_rmse)
# Convert to a DataFrame and find the optimal alpha
alpha_fit = pd.DataFrame({'alpha': alphas, 'RMSE': ses_rmse_values})
alpha_min = alpha_fit.loc[alpha_fit['RMSE'].idxmin()]
# Plot RMSE vs. alpha
plt.figure(figsize=(10, 6))
plt.plot(alpha_fit['alpha'], alpha_fit['RMSE'], marker='o')
plt.scatter(alpha_min['alpha'], alpha_min['RMSE'], color='blue', s=100)
plt.xlabel('Alpha')
plt.ylabel('RMSE')
plt.title('RMSE vs. Alpha for SES')
plt.show()
The plot above shows that alpha should be 0.1. Thus, we will fit the SES model using the training set and setting the alpha parameter to 0.1. The following code performs this task and generates the model forecast graph below.
# Refit the optimal model with the best alpha
best_alpha = alpha_min['alpha']
ses_model = SimpleExpSmoothing(ses_train).fit(smoothing_level=best_alpha, optimized=False)
ses_forecast = ses_model.forecast(len(ses_test))
# Plotting results
plt.figure(figsize=(10, 6))
plt.plot(ses_train, label='Train')
plt.plot(ses_test, label='Test')
plt.plot(ses_model.fittedvalues, label='Fitted')
plt.plot(ses_test.index, ses_forecast, label='Forecast', color='red')
plt.xlabel('Time')
plt.ylabel('Log Sales')
#plt.title('Simple Exponential Smoothing')
plt.legend()
plt.show()
The model captures the overall trend in the data. However, it fails to capture the sudden changes around the average performance. Predictions for the test set are flat.
Forecast performance evaluation using SES¶
To evaluate the forecast performance based on the SES method, we compute the RMSE metric. See below.
# Performance evaluation
ses_train_rmse = np.sqrt(mean_squared_error(ses_train, ses_model.fittedvalues))
ses_test_rmse = np.sqrt(mean_squared_error(ses_test, ses_forecast))
print(f'Train RMSE: {ses_train_rmse}')
print(f'Test RMSE: {ses_test_rmse}')
Train RMSE: 0.3603028868685023 Test RMSE: 0.39489577005480325
The output above shows that the RMSE for our test set is 0.392, while the RMSE for the training set is 0.358. Comparing the RMSE values from the training and test sets suggests that there is no overfitting issue.
Sales forecast using SES¶
Finally, we will produce the sales forecasts for the next 12 months. The following code chunk generates the forecasts from SES. Since the model uses the logged sales, we take the antilogarithm using the np.exp() function in Python to obtain the final sales forecasts.
# Estimate the SES model using the entire dataset and generate forecasts for the next 12 months
ses_final_model = SimpleExpSmoothing(data['LSales']).fit(smoothing_level=best_alpha, optimized=False)
ses_final_forecast = ses_final_model.forecast(12)
# Calculate Mean Squared Error for the entire model
ses_residuals = ses_final_model.resid
ses_mse = np.mean(np.square(ses_residuals))
# Calculate prediction intervals
forecast_steps = 12
alpha = 0.05 # 95% prediction intervals
ses_forecast_se = np.sqrt(ses_mse) * np.sqrt(np.arange(1, forecast_steps + 1)) # Forecast standard error
ses_ci_lower = ses_final_forecast - 1.96 * ses_forecast_se
ses_ci_upper = ses_final_forecast + 1.96 * ses_forecast_se
# Convert forecasts and intervals to original scale if log transformation was used
ses_final_forecast_sales = np.exp(ses_final_forecast)
ses_ci_lower_sales = np.exp(ses_ci_lower)
ses_ci_upper_sales = np.exp(ses_ci_upper)
# Create a DataFrame for the forecasted sales with prediction intervals
forecast_df = pd.DataFrame({
'Forecast': ses_final_forecast_sales,
'Lower Bound': ses_ci_lower_sales,
'Upper Bound': ses_ci_upper_sales
})
# Display the DataFrame
print('SES model sales forecasts for the next 12 months:')
print(forecast_df)
SES model sales forecasts for the next 12 months:
Forecast Lower Bound Upper Bound
122 201.651411 97.764511 415.931006
123 201.651411 72.434287 561.381812
124 201.651411 57.545292 706.631080
125 201.651411 47.398129 857.909207
126 201.651411 39.951832 1017.807929
127 201.651411 34.232027 1187.872718
128 201.651411 29.697785 1369.236528
129 201.651411 26.018791 1562.843215
130 201.651411 22.979531 1769.543979
131 201.651411 20.432327 1990.144945
132 201.651411 18.272080 2225.433070
133 201.651411 16.421708 2476.191368
The output above shows the point forecasts for the next 12 months alongside the prediction intervals. The first column shows the point forecast , and the second and third columns show the lower and upper bounds of the prediction interval at a 0.95 probability. For example, the first row suggests that predicted sales for the next month are 202.731. With a 0.95 probability, the next month’s sales will fall within the range of [98.738, 416.252].
These are sales forecasts from the SES method. Now, we proceed to the second method: Holt–Winters.
Holt-Winters method¶
Similar to SES model building, in the Holt–Winters method, we first create a training set, as shown in the following code chunk.
from statsmodels.tsa.holtwinters import ExponentialSmoothing
hw_train = data['LSales'][:96]
Model fitting with the Holt-Winters method¶
Before we predict future values, we will need to fit the Holt–Winters model to the data. If we are dealing with many time-series data simultaneously, we may simply want to call the Holt–Winters function in Python and let the software figure out the tuning parameters of the model on its own. We also have the opportunity to tune the fit manually:
alpha: The parameter for the level equation. A higher alpha puts more weight on the most recent observations.
beta: The parameter for the trend equation. A higher beta means the trend slope is more dependent on recent trend slopes.
gamma: The parameter for the seasonality equation. A higher gamma puts more weight on the most recent seasonal observations.
For further information about model parameters, please see the corresponding chapter of the book.
For model fitting, we may have two approaches:
HW1: This refers to the model where Python automatically finds the best alpha, beta, and gamma values.
HW2: This refers to the model where we specify alpha, beta, and gamma to be 0.2, 0.1, and 0.2, respectively. We plot the fitted values versus the raw data to see the quality of the model fit. We may want to try different values for alpha, beta, and gamma to see how the fit changes.
The below code obtains the model fit using both approaches, HW1 and HW2.
warnings.filterwarnings("ignore")
# Fitting Holt-Winters model with default smoothing parameters
HW1 = ExponentialSmoothing(hw_train, seasonal='add', seasonal_periods=12).fit()
# Custom Holt-Winters' fitting with specified smoothing parameters
HW2 = ExponentialSmoothing(hw_train, seasonal='add', seasonal_periods=12).fit(smoothing_level=0.2, smoothing_slope=0.1, smoothing_seasonal=0.1)
# Visualizing the fits
plt.figure(figsize=(10, 6))
plt.plot(hw_train, label='Observed',color='black')
plt.plot(HW1.fittedvalues, label='HW1 Fitted', linestyle='--', color='blue')
plt.plot(HW2.fittedvalues, label='HW2 Fitted', linestyle='--', color='red')
plt.xlabel('Time')
plt.ylabel('LnSales')
plt.legend()
plt.show()
The solid black line in this plot is the log-sales data. The blue dashed line shows the model fit based on HW1, while the red dashed line represents the model fit based on HW2. How well does the Holt–Winters method capture the past patterns in the data? The actual sales data exhibit frequent fluctuations. Both HW1 and HW2 capture the overall level and trend of the sales data reasonably well but do not fully capture the high variability in the data.
Next, we will evaluate forecast performance using the Holt-Winters method with the HW1 model fit.
Forecast performance evaluation using the Holt-Winters method¶
First, we generate the predicted sales in the test set based on the HW1 model fit. To this end, we use the following code and obtain the forecast plot below.
# Holt-Winters Model Forecasts with Prediction Intervals
forecast_steps = 26
HW1_pred = HW1.forecast(steps=forecast_steps)
# Calculate Mean Squared Error
residuals = HW1.resid
mse = np.mean(np.square(residuals))
# Calculate prediction intervals
alpha = 0.05 # 95% prediction intervals
forecast_se = np.sqrt(mse) * np.sqrt(np.arange(1, forecast_steps + 1)) # Forecast standard error
ci_lower = HW1_pred - 1.96 * forecast_se
ci_upper = HW1_pred + 1.96 * forecast_se
# Visualizing the forecast with prediction intervals
plt.figure(figsize=(10, 6))
plt.plot(hw_train, label='Observed', color='black')
plt.plot(HW1.fittedvalues, label='Fitted Sales', linestyle='--', color='blue')
plt.plot(HW1_pred, label='Forecast', color='red')
plt.fill_between(HW1_pred.index, ci_lower, ci_upper, color='pink', alpha=0.3, label='Prediction Interval')
#plt.title('Holt-Winters Forecast with Prediction Intervals')
plt.xlabel('Time')
plt.ylabel('LnSales')
plt.legend()
plt.show()
The solid black line is the observed log-sales data. The blue dashed line is the estimated sales based on the HW1 model fit. The red solid line represents the predicted sales over the next quarter, along with the upper and lower bounds of the prediction interval with 0.95 probability. As the plot suggests, forecasts mimic the average pattern in the data, and the prediction intervals spread wildly outward.
Multiplicative seasonality
When we use the Holt–Winters method to fit the data, we may also want to tune the behaviour of the seasonality component. The standard Holt–Winters method uses additive seasonality, which assumes that the amplitude of any seasonality components is relatively constant throughout the series. However, if we use multiplicative seasonality, we allow the seasonal variations to grow with the overall level of the data. To see how that works, we will fit the model, generate forecasts, and compare the results to our additive fit of HW1. The following code produces the forecast plot based on the choice of multiplicative seasonality in the Holt–Winters method. As the plot suggests, the confidence intervals spread wildly outward. For this dataset, multiplicative seasonality fitting does not appear to be the way to go.
# Holt-Winters' Model with Multiplicative Seasonality
# Fit Holt-Winters model with multiplicative seasonality
HW3= ExponentialSmoothing(hw_train, seasonal='mul', seasonal_periods=12).fit()
# Summary of the model
print(HW3.summary())
# Forecasting with prediction intervals
forecast_steps = 26
HW3_forecast = HW3.forecast(steps=forecast_steps)
# Calculate Mean Squared Error
residuals_hw3 = HW3.resid
mse_hw3 = np.mean(np.square(residuals_hw3))
# Calculate prediction intervals
alpha = 0.05 # 95% prediction intervals
forecast_se_hw3 = np.sqrt(mse_hw3) * np.sqrt(np.arange(1, forecast_steps + 1)) # Forecast standard error
ci_lower_hw3 = HW3_forecast - 1.96 * forecast_se_hw3
ci_upper_hw3 = HW3_forecast + 1.96 * forecast_se_hw3
# Visualizing the forecast with prediction intervals
plt.figure(figsize=(10, 6))
plt.plot(hw_train, label='Observed', color='black')
plt.plot(HW3.fittedvalues, label='Fitted Sales', linestyle='--', color='blue')
plt.plot(HW3_forecast, label='Forecast', color='red')
plt.fill_between(HW3_forecast.index, ci_lower_hw3, ci_upper_hw3, color='pink', alpha=0.3, label='Prediction Interval')
#plt.title('Holt-Winters Forecast with Multiplicative Seasonality')
plt.xlabel('Time')
plt.ylabel('LnSales')
plt.legend()
plt.show()
ExponentialSmoothing Model Results
================================================================================
Dep. Variable: LSales No. Observations: 96
Model: ExponentialSmoothing SSE 10.313
Optimized: True AIC -186.167
Trend: None BIC -150.266
Seasonal: Multiplicative AICC -179.281
Seasonal Periods: 12 Date: Wed, 22 May 2024
Box-Cox: False Time: 20:49:09
Box-Cox Coeff.: None
=================================================================================
coeff code optimized
---------------------------------------------------------------------------------
smoothing_level 1.4901e-08 alpha True
smoothing_seasonal 0.000000 gamma True
initial_level 5.0799331 l.0 True
initial_seasons.0 1.0419272 s.0 True
initial_seasons.1 1.0727155 s.1 True
initial_seasons.2 1.0626704 s.2 True
initial_seasons.3 1.0563466 s.3 True
initial_seasons.4 1.0475162 s.4 True
initial_seasons.5 1.0445707 s.5 True
initial_seasons.6 1.0460946 s.6 True
initial_seasons.7 1.0318329 s.7 True
initial_seasons.8 1.0742244 s.8 True
initial_seasons.9 1.0245300 s.9 True
initial_seasons.10 1.0244080 s.10 True
initial_seasons.11 1.0603786 s.11 True
---------------------------------------------------------------------------------
For this dataset, multiplicative seasonality fitting provides similar forecasts to the HW1 approach. Therefore, we turn to the HW1 approach which is less sophisticated.
Having conducted these explorations, we calculate the RMSE of predicted sales using HW1 approach for both the training set and the test set, as shown in the following code chunk.
# RMSE of HW1
forecast_steps = 26
HW1_pred = HW1.forecast(steps=forecast_steps)
# Calculate Mean Squared Error
residuals = HW1.resid
mse = np.mean(np.square(residuals))
# Calculate Mean Squared Error for the training set
train_residuals = HW1.resid
train_mse = np.mean(np.square(train_residuals))
# Calculate RMSE for the training set
train_rmse = np.sqrt(train_mse)
# Calculate RMSE for the test set (observations 97 to 122)
test_set = data['LSales'][96:96 + forecast_steps]
test_rmse = np.sqrt(np.mean(np.square(test_set - HW1_pred)))
# Print RMSE for both training and test sets
print(f"RMSE of HW1 on the training set: {train_rmse}")
print(f"RMSE of HW1 on the test set: {test_rmse}")
RMSE of HW1 on the training set: 0.3277673603047565 RMSE of HW1 on the test set: 0.40656084859722563
The model’s output above indicates that the RMSE of HW1 on the training set is 0.326, and the RMSE of HW1 on the test set is 0.403. It is worth noting that if the forecast accuracy metrics from the training and test sets are close to each other, the model is considered robust in terms of fit. Next, we will use the Holt–Winters method with the HW1 approach to generate sales forecasts.
Sales forecast using Holt-Winters method¶
We will generate sales forecasts based on the Holt-Winters with the HW1 approach including prediction intervals for the next 12 months. Note that these are the forecasts for the sales data in logarithm. Therefore, we need to convert the logged sales forecasts to sales forecasts in level. We do this by taking the exponentiation of the forecasts. See below.
HW_final = ExponentialSmoothing(data['LSales'], seasonal='add', seasonal_periods=12).fit()
forecast_steps = 12
HW_final_forecast = HW_final.forecast(steps=forecast_steps)
# Calculate Mean Squared Error
residuals = HW_final.resid
mse = np.mean(np.square(residuals))
# Calculate prediction intervals
alpha = 0.05 # 95% prediction intervals
forecast_se = np.sqrt(mse) * np.sqrt(np.arange(1, forecast_steps + 1)) # Forecast standard error
ci_lower = HW_final_forecast - 1.96 * forecast_se
ci_upper = HW_final_forecast + 1.96 * forecast_se
# Convert forecasts and intervals to original scale if log transformation was used
HW_final_forecast_sales = np.exp(HW_final_forecast)
ci_lower_sales = np.exp(ci_lower)
ci_upper_sales = np.exp(ci_upper)
# Create a data frame for the forecasts and prediction intervals
hw_forecasts_df= pd.DataFrame({
'Forecast': HW_final_forecast_sales,
'Lower Bound': ci_lower_sales,
'Upper Bound': ci_upper_sales
})
# Display the DataFrame
hw_forecasts_df
# Visualize the HW forecast with prediction intervals
#plt.figure(figsize=(10, 6))
#plt.plot(data['Sales'], label='Observed', color='black')
#plt.plot(HW_final.fittedvalues, label='Fitted Sales', linestyle='--', color='blue')
#plt.plot(HW_final_forecast_sales, label='Forecast', color='red')
#plt.fill_between(HW_final_forecast_sales.index, ci_lower_sales, ci_upper_sales, color='pink', alpha=0.3, label='Prediction Interval')
#plt.title('Holt-Winters Forecast with Prediction Intervals')
#plt.xlabel('Date')
#plt.ylabel('Sales')
#plt.legend()
#plt.show()
| Forecast | Lower Bound | Upper Bound | |
|---|---|---|---|
| 122 | 217.086524 | 110.818376 | 425.259425 |
| 123 | 214.053783 | 82.707088 | 553.991479 |
| 124 | 202.468777 | 63.177539 | 648.863602 |
| 125 | 194.612795 | 50.714158 | 746.815901 |
| 126 | 204.248145 | 45.412889 | 918.622574 |
| 127 | 183.765893 | 35.396658 | 954.042138 |
| 128 | 208.094336 | 35.127326 | 1232.751169 |
| 129 | 192.703815 | 28.769344 | 1290.775355 |
| 130 | 193.050313 | 25.680724 | 1451.221685 |
| 131 | 215.327597 | 25.683130 | 1805.308556 |
| 132 | 205.573034 | 22.102519 | 1912.011541 |
| 133 | 227.785387 | 22.178652 | 2339.465076 |
The output above shows the point forecasts for the next 12 months alongside the prediction intervals. The first column shows the point forecast , and the second and third columns show the lower and upper bounds of the prediction interval at a 0.95 probability. For example, the first row suggests that predicted sales for the next month are 218.152. With a 0.95 probability, the next month’s sales will fall within the range of [111.837, 425.533].
Next, we will implement the ARIMA modelling approach.
ARIMA model¶
We will build an ARIMA model following the procedure summarized below:
Perform unit root tests to check for non-stationarity in the data and perform differencing if necessary.
Obtain the ACF and PACF plots to determine the order of lags and hence the specification of ARIMA.
Split the data into training and test sets, and estimate an ARIMA model using the training set and predict sales using the test set. We will evaluate forecast accuracy via the RMSE metric.
Generate the sales forecast for the next twelve months.
Similar to the previous models, we will use the log-transformed sales data. Our first step is to conduct the unit root tests to exanine the stationarity of the data. The visual inspection of the sales performance plot suggests that there is some variation in the data. To reduce the variation, we will take the logarithm:
Stationarity test
As introduced in the chapter, there are multiple tests for testing stationarity, including the ADF, KPSS, and Phillips-Perron test. In this application, we demonstrate the procedure with the ADF test. Recall that under the ADF test, the null and alternative hypotheses are:
- $H_{0}$: The data is not stationary.
- $H_{a}$: The data is stationary.
If the p-value from the ADF test is smaller than a certain significance level (e.g .05), we should reject the null hypothesis and conclude that the variable is stationary.
To run the ADF test in Python, we use the adfuller function.
# Augmented Dickey-Fuller test on Log-transformed Sales
result = adfuller(data['LSales'].dropna())
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
ADF Statistic: -7.903300 p-value: 0.000000
We reject the null hypothesis because the obtained p-value, 0, is smaller than the significance level, .05. Therefore, the ADF test result suggests that the series is stationary.
ACF and PACF analysis
A systematic way to determine the order of lags for the autoregressive (AR) and moving average (MA) components of ARIMA model is to plot and inspect the ACF and PACF graphs. The following code chunk generates (i) the ACF plot and (ii) the PACF plot, simultaneously:
# Time series display for Sales (Trend, ACF, PACF)
plt.figure(figsize=(10, 6))
plt.subplot(211)
plot_acf(data['LSales'].dropna(), ax=plt.gca())
plt.subplot(212)
plot_pacf(data['LSales'].dropna(), ax=plt.gca())
plt.show()
In the ACF and PACF plots, the dashed horizontal lines represent the upper and lower bounds of the critical region (i.e. 95% confidence level). The lag order of the AR and MA components is identified by the behaviour of the PACF and ACF plots. We find that both ACF and PACF have a cut-off at lag 1. This suggests that we should probably take a lag order of 1 for both MA and AR components. Given that we did not take the difference of the sales series, the specification of the non-seasonal part of our model is ARIMA (1,0,1).
Furthermore, the sales data show a recurring pattern and the ACF plot displays a sinusoidal waves. This kind of pattern suggests the existence of seasonality in the data. Therefore, the best specification would be seasonal ARIMA (1,0,1) (1,0,1). Note that at this point we assume the same number of lags for the AR and MA components in the seasonal part. However, we will explore this specification in more detail below when we discuss the model estimation.
Data split¶
To examine the ARIMA model’s predictive ability, we need to split the data into training (in-sample) and test (out-of-sample) sets.
To do this, we apply the 80/20 rule, using the first 80% of the observations for the training set and the remaining 20% for the test set. Given that we have 121 observations in total, we should use the first 96 observations for the training set, and the remaining 25 observations for the test set. The following code chunk shows how we implement the data split in RPython.
# Preparing training and testing sets for ARIMA modeling
arima_train = data['LSales'][:96]
arima_test = data['LSales'][96:121]
Model estimation¶
Recall that the ACF and PACF analysis suggests the seasonal ARIMA (1,0,1) (1,01) model specification. To estimate this model using the training set, we use the Python function ARIMA:
import warnings
warnings.filterwarnings("ignore")
# Estimating the ARIMA model using the training set
model = ARIMA(arima_train, order=(1, 0, 1), seasonal_order=(1, 0, 1, 12))
fit_arima = model.fit()
print(fit_arima.summary())
SARIMAX Results
========================================================================================
Dep. Variable: LSales No. Observations: 96
Model: ARIMA(1, 0, 1)x(1, 0, 1, 12) Log Likelihood -22.220
Date: Wed, 22 May 2024 AIC 56.440
Time: 20:51:36 BIC 71.826
Sample: 0 HQIC 62.659
- 96
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 5.3362 0.059 90.300 0.000 5.220 5.452
ar.L1 0.2578 0.267 0.965 0.334 -0.266 0.781
ma.L1 0.1844 0.265 0.696 0.487 -0.335 0.704
ar.S.L12 -0.8237 0.424 -1.943 0.052 -1.655 0.007
ma.S.L12 0.9989 28.516 0.035 0.972 -54.891 56.889
sigma2 0.0865 2.435 0.036 0.972 -4.686 4.859
===================================================================================
Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 11.90
Prob(Q): 0.96 Prob(JB): 0.00
Heteroskedasticity (H): 1.18 Skew: -0.16
Prob(H) (two-sided): 0.64 Kurtosis: 4.69
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
This output shows:
the estimated model parameters p, q , P and Q (displayed as 'ar.L1', 'ma.L1', 'ar.S.L12' and 'ma.S.L12', respectively) and their standard errors. For example, p is estimated to be 0.2578 and its standard error is 0.267.
the residuals variance $\sigma^2$ (displayed as 'sigma2').
the model fit statistics (log likelihood, AIC, AICc and BIC).
the model diagnostics such as Ljung-Box and Jarque-Bera tests.
Although we relied on the ACF and PACF analysis for model identification, we may want to estimate ARIMA models with different specifications and compare model performances based on model fit statistics such as AIC and BIC.
There are also other Python functions that can automatically pick the specifications with the lowest AIC and BIC--for example, auto_arima function from the pmdarima. However, we should keep in mind that model identification with an automated procedure may sometimes provide misleading results. For instance, some firms might operate within a certain cycle and would want to evaluate sales using a specific order of lags.
Forecast performance evaluation of seasonal ARIMA¶
Next, we evaluate the forecast performance of the model using the fitted model parameters.
import seaborn as sns
# Predict sales based on seasonal ARIMA model estimations
predict_arima = fit_arima.get_forecast(steps=25)
predict_arima_pred = predict_arima.predicted_mean
# Plot the actual and predicted sales on the same graph using seaborn
data_plot = pd.DataFrame({
'Month': arima_test.index,
'Actual Sales': arima_test,
'Predicted Sales': predict_arima_pred
})
# Melt the DataFrame to long format for seaborn
data_plot_long = pd.melt(data_plot, id_vars=['Month'], var_name='Type', value_name='Sales')
# Plot using seaborn
plt.figure(figsize=(10, 6))
sns.lineplot(data=data_plot_long, x='Month', y='Sales', hue='Type', style='Type', dashes=[(2, 2), (1, 0)])
plt.show()
This plot suggests that the predicted sales based on the seasonal ARIMA model can mimic (to a certain extent) the general pattern of actual sales in the test set. To determine which method does a relatively more accurate job, we will calculate and compare the root mean squared error (RMSE) of both predictions.
Sales forecast using seasonal ARIMA¶
The below code computes the RMSE metric for the seasonal ARIMA model.
# Calculate RMSE for the training set
train_pred_arima = fit_arima.fittedvalues
rmse_train_arima = mean_squared_error(arima_train, train_pred_arima, squared=False)
# Calculate RMSE for the test set
rmse_test_arima = mean_squared_error(arima_test, predict_arima_pred, squared=False)
# Print RMSE for both training and test sets
print(f"RMSE of ARIMA on the training set: {rmse_train_arima}")
print(f"RMSE of ARIMA on the test set: {rmse_test_arima}")
RMSE of ARIMA on the training set: 0.3053227567058221 RMSE of ARIMA on the test set: 0.3525117588639769
The RMSEs of our seasonal ARIMA on the training set and test set are 0.305 and 0.353, respectively.
Finally, we will re-estimate the seasonal ARIMA model using the full data set, and then generate twelve-month sales forecasts.
# Estimating the ARIMA model using the entire dataset and generating forecasts for the next 12 months
fit_arima_final = ARIMA(data['LSales'], order=(1, 0, 1), seasonal_order=(1, 0, 1, 12)).fit()
forecast_arima_final = fit_arima_final.get_forecast(steps=12)
forecast_arima_final_mean = forecast_arima_final.predicted_mean
# Print final forecasts
print(forecast_arima_final_mean)
122 5.150800 123 5.242018 124 5.359608 125 5.369407 126 5.299007 127 5.406623 128 5.320422 129 5.370669 130 5.407167 131 5.289700 132 5.285022 133 5.197894 Name: predicted_mean, dtype: float64
The output above shows the point forecasts for the next twelve months. We can also plot these forecasts as follows:
# Plotting sales forecasts using SARIMA
plt.plot(forecast_arima_final_mean, label='Forecast')
plt.legend()
plt.show()
Let's put the above one into the historical sales graph to extend it into future.
import matplotlib.pyplot as plt
import pandas as pd
forecast_steps = 12
# Convert the historical index to a pandas datetime index if it's not already one
data['LSales'].index = pd.to_datetime(data['LSales'].index)
# Plot historical sales
plt.figure(figsize=(12, 7))
plt.plot(data['LSales'], label='Historical Sales', color='skyblue')
# The last date of historical data
last_date = data['LSales'].index[-1]
# The frequency of the historical data is used to create a correct forecast index
if isinstance(data['LSales'].index, pd.DatetimeIndex):
# Get the frequency of the historical data
freq = pd.infer_freq(data['LSales'].index)
# Generate forecast index starting after the last historical date
forecast_index = pd.date_range(start=last_date, periods=forecast_steps + 1, freq=freq)[1:]
else:
# If the index is not a datetime index, generate a range of integers for plotting
forecast_index = range(len(data['LSales']), len(data['LSales']) + forecast_steps)
# Add forecasted sales to the plot
plt.plot(forecast_index, forecast_arima_final_mean, label='Forecast', color='blue')
# Get confidence intervals for the forecast
confidence_intervals = fit_arima.get_forecast(steps=forecast_steps).conf_int()
# Plot the confidence intervals
plt.fill_between(forecast_index, confidence_intervals.iloc[:, 0], confidence_intervals.iloc[:, 1], color='blue', alpha=0.1)
plt.legend()
plt.show()
Note that these are the forecasts for the sales data in logarithm. Therefore, we need to convert the logged sales forecasts to sales forecasts in level. We do this by taking the exponentiation of the forecasts. Also, we need to compute the prediction intervals for these forecasts. See below.
# Convert final forecasts to original scale
forecast_sales_arima_final = np.exp(forecast_arima_final_mean)
# Get confidence intervals for the forecast
confidence_intervals = forecast_arima_final.conf_int()
# Convert prediction intervals to original scale
ci_lower_sales = np.exp(confidence_intervals.iloc[:, 0])
ci_upper_sales = np.exp(confidence_intervals.iloc[:, 1])
# Create a DataFrame for the forecasted sales with prediction intervals
forecast_df = pd.DataFrame({
'Forecast': forecast_sales_arima_final,
'Lower Bound': ci_lower_sales,
'Upper Bound': ci_upper_sales
})
# Display the DataFrame
print(forecast_df)
Forecast Lower Bound Upper Bound 122 172.569446 92.423860 322.213482 123 189.051198 98.284487 363.641879 124 212.641659 110.438339 409.427335 125 214.735478 111.521960 413.473055 126 200.137979 103.940688 385.366032 127 222.877559 115.750374 429.151152 128 204.470058 106.190528 393.707474 129 215.006616 111.662648 413.995600 130 222.998900 115.813698 429.383660 131 198.283950 102.985715 381.766780 132 197.358510 102.583116 379.695832 133 180.890817 94.023530 348.013821
These are our final sales forecasts using the seasonal ARIMA modelling approach.
So far, we have executed three different models (SES, Holt–Winters, and ARIMA) and obtained three different sets of forecasts. How should we decide which model to choose?
Model comparison¶
We compare the forecast performance of these models based on the RMSE metric from the test set. . The following code summarizes the RMSE output and displays it in a table.
#RMSE metrics from the above
rmse_ses = 0.392
rmse_holtwinters = 0.403
rmse_arima = 0.353
# Create a DataFrame to summarize the RMSE values
table_rmse_summary = pd.DataFrame({
'Method': ['SES', 'Holt-Winters','ARIMA'],
'RMSE': [ rmse_ses, rmse_holtwinters,rmse_arima]
})
# Display the DataFrame
print(table_rmse_summary)
Method RMSE 0 SES 0.392 1 Holt-Winters 0.403 2 ARIMA 0.353
Which of these models should Silverline choose, and why?
The RMSE values suggest that the seasonal ARIMA model outperforms the other two models. Recall the final forecasts from the seasonal ARIMA model (Figure 9.30).
All in all, using the seasonal ARIMA model’s forecasts for this store, Silverline can try to predict the demand for the other store struggling with the stockout problems. Further investigation may be warranted, as there are doubtless other sophisticated models which include further significant drivers of store demand (e.g. price, promotion) that would improve these forecasts.
Summary¶
In this tutorial, we demonstrated the application of the demand forecasting models (SES, Holt-Winters, and ARIMA) in Python to a kitchen appliance company from an emerging market. More specifically, we emphasized the importance of model validation and showed that once a forecasting model is developed to fit past data well, it can be used to predict what may happen in the future.
The modelswe covered in this tutorial are considered less sophisticated than some of the advanced machine learning models (e.g. regression trees, neural networks, support vector machines), but their performance sometimes may be surprisingly better (Makridakis et al., 2018). Therefore, our recommendation is that their predictive performance should always be compared with that of more complex models available to the user . It is important to note that in using all these predictive models we make an assumption that consumers’ past purchases are good indicators of future demand and that changes in purchasing context (e.g. shopping companions, consumer’s mood) do not alter consumer preferences. However, research in consumer psychology often shows that in some situations consumers are inconsistent in their preferences (Israeli and Avery, 2017). For example, consumers sometimes may seek variety in their purchases and thus may not show consistent buying behaviour. In the same vein, consumers may make impulse purchases rather than planned purchases. This kind of behaviour may be difficult to predict based on past data. Moreover, social context (e.g. shopping in the presence of friends) may influence consumers’ buying habits. Therefore, it is recommended that marketing analysts understand the conceptual underpinnings of consumer purchase behaviour prior to embarking on a demand forecasting task, and investigate the past data of the product accordingly.